Antiferromagnetism in hexagonal graphene structures: rings vs dots 



o 

(N 



M. Grujic/'B M. Tadic/'Q and F. M. Peeters^-I 

^School of Electrical Engineering, University of Belgrade, P.O. Box 3554, IUZO Belgrade, Serbia 
'^Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020 Antwerp, Belgium 

The mean-field Hubbard model is used to investigate the formation of the antiferromagnetic 
phase in hexagonal graphene rings with inner zigzag edges. The outer edge of the ring was taken 
to be either zigzag or armchair, and we found that both types of structures can have a larger 
antiferromagnetic interaction as compared with hexagonal dots. This difference could be partially 
ascribed to the larger number of zigzag edges per unit area in rings than in dots. Furthermore, edge 
states localized on the inner ring edge are found to hybridize differently than the edge states of dots, 
which results in important differences in the magnetism of graphene rings and dots. The largest 
staggered magnetization is found when the outer edge has a zigzag shape. However, narrow rings 
with armchair outer edge are found to have larger staggered magnetization than zigzag hexagons. 
The edge defects are shown to have the least effect on magnetization when the outer ring edge is 
armchair shaped. 
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Recent tremendous progress in graphene research is 
driven by its remarkable properties, e.g. high crystalline 
quality, high electron mobility, lack of a band gap, and 
a minimal possible thickness, to name a fewii The men- 
tioned properties are advantageous for various applica- 
tions of graphene, such as piezoelectric devices;^ superca- 
pacitorSf^ photodetectorsyi and field-effect transistors<^ 
Furthermore, it has been predicted that graphene struc- 
tures could exhibit magnetic ordering which is potentially 
advantageous for spintronic applications This effect is 
essentially related to either the global or local imbalance 
of sublattice atoms in bipartite lattices. An imbalance 
might give rise to zero energy states in the electron spec- 
trum. These states are localized near the zigzag edges or 
vacancies, and along with the repulsive electron-electron 
(e-e) interaction could eventually lead to a spin polariza- 
tion of the ground state of the systemji Furthermore, the 
spins on the same sublattice are found to exhibit ferro- 
magnetic coupling along the graphene edges, whereas the 
spins on different sublattices along the graphene edges 
couple antiferromagnetically. 

In theory, magnetic ordering has been demonstrated 
for graphene flakeS)^ nanoribbonSfiS and vacancies in 
bulk graphenciii On the other hand, experimental re- 
ports on magnetism in graphene structures are rare and 
conflicting. They range from the detection of ferromag- 
netic or antiferromagnetic orderingi^^— to measurements 
of defect-induced paramagnetism! ^^'^^ Magnetic order- 
ing was even found to be preserved at room tempera- 
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The essential cause of magnetism in graphene 
is the existence of a peak in the density of nonbond- 
ing edge states near the Fermi energy. However, due 
to the high reactivity of these states, magnetism might 
be strongly suppressed?^ Several theoretical studies of- 
fered explanations for a diversity of phenomena related 
to magnetic ordering and its suppression, which might 
occur by means of nonmagnetic edge passivation, edge 
reconstruction, or vanishing of spin correlations with in- 



creasing temperaturcj Hence, in order to experimen- 
tally detect magnetic ordering graphene samples should 
be kept under rigorously controlled conditions. Yet, var- 
ious applications of this effect have been proposed. They 
involve half-metallicity with electrically controlled spin 
propagation;^ defect induced spin filtering^ and spin 
logic devices.— 1^ 

In this report we employ the mean-field Hubbard 
model to study the formation of local magnetic moments 
in hexagonal graphene rings. Our aim is to explore how 
magnetic ordering is affected by the ring size and the 
edge type. In order to identify different hexagonal rings, 
we introduce the following notation which might be vi- 
sually aiding. We assume that the type of the inner ring 
edge is zigzag, and that N unit hexagons are adjacent 
to this boundary. The outer ring edge is assumed to be 
comprised of either M diniers if it is of armchair type, or 
M unit hexagons if it is of zigzag type. Therefore, the 
ring is denoted by Af : N . As an example, consider the 
ring shown in Fig. 1, which is assumed to be formed out 
of the hexagonal dot with armchair edge, which contains 
seven dimers at each side of the hexagon, as shown in 
Fig. 1(a). The ring is formed when the carbon atoms 
around the center of the dot are removed, as depicted 
in Figs. 1(a) and (b). Potentially these exotic struc- 
tures could be manufactured via substitutional doping of 
boron nitride nanostructures with carbon.— Because the 
edge of the removed dot has four unit hexagons at each 
side, the ring is denoted as 7ac ■ ^zg- The distribu- 
tions of the magnetic moments in the graphene rings will 
be compared with the magnetic moment distributions in 
the hexagonal graphene dots. Those dots are assumed to 
have zigzag edges, and are labeled by Nzgi where TV has 
the same meaning as the symbol M for the rings. 

Magnetic ordering of a graphene structure is governed 
by Lieb's theorem.— It states that the total ground-state 
spin of a bipartite lattice with repulsive e-e interaction 
as described by the Hubbard model equals half of the 
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difference of the sublattice sites. For symmetrical struc- 
tures, this rule is related to the arrangement of the carbon 
atoms with respect to lines of reflection symmetry in the 
graphene plane: if there is a symmetry line which does 
not intersect any of the carbon atoms the total ground 
state spin is zero; otherwise there exists a finite magnetic 
moment. All the hexagonal rings analyzed here possess 
such a symmetry, thus their total magnetization equals 
zero, unlike triangular rings which display a ferrimag- 
netic phase<^ However, Lieb's theorem does not dictate 
the distribution of the local magnetic moments or the 
lack of zero energy states. Furthermore, the number of 
zero-energy states in the analyzed ring is an integer mul- 
tiple of six, which is a consequence of the Cqv symmetry 
of the ring. In the 9ac '■ IOzg rhig, six zero energy states 
are found, which agrees with graph theory, and which is a 
topological property related to the nonperfect matching 
of the pz orbitalsii 

The Hubbard Hamiltonian 

H = Ho + Hi, (1) 

is employed to compute the distribution of magnetic mo- 
ments. Hq is the noninteracting part, which represents 
the nearest neighbor tight-binding Hamiltonian, and is 
given by 

Ho^-t J2 4c,.: (2) 

<i,j>,(J 

where Cja and Cj^ are the annihilation and creation op- 
erators, respectively, and t denotes the hopping integral. 
The interacting part Hi describes the e-e interaction 

Hi (n.^t(?^•4) + nii{nif) - (ni^){nii)) , (3) 

i 

where rii^ = cl^Ci„ is the number operator, and U de- 
notes the on-site Coulomb repulsion energy for each pair 
of electrons with the opposite spins orbiting the same 
atom^i Equation ([3]) is obtained within the mean-field 
approximation, which assumes that the spin-up (spin- 
down) electrons interact with the average density of spin- 
down (spin-up) electrons on a particular atomic site. 

In our calculations, we take t = 2.7 eV and U = 1.2t^ 
We note that there is no consensus on the actual value of 
the strength of the Coulomb interaction to be used in the 
Hubbard model in graphene. Recent density functional 
theory (DFT) calculations came up with a value closer to 
U = 3.4^1^ However, having in mind that the mean-field 
approximation can overestimate the tendency for mag- 
netic order for large U^^ we chose the more conservative 
value of [/ — 1.2t. The solution is then obtained by 
means of a self-consistent procedure which starts from 
an initial distribution of the spins, and ends when the 
maximum change of the electron density over the atomic 
sites drops below 10~^. When the self-consistent spin 
densities are determined, the magnetic moment per site 
nii is computed as 

= «') = (Kt>-("4))/2- (4) 




FIG. 1: (Color online) (a) The dot (red color) with four atoms 
at the zigzag edge removed from the larger dot (black color) 
which has seven dimers at the dot edge, (b) The formed ring 
is labeled by 7ac '■ 4zg. (c) The distribution of magnetic 
moments in the 9ac '■ Nza ring shown in a sextant of the ring 
for A'' taking the values 3, 4, 5, 9, 10 and 11. The majority 
spin is labeled by both orientation and color of a triangle 
centered at an atomic site. The local magnetic moment value 
is proportional to the color intensity. 

For the antiferromagnetic order parameter we take the 
staggered magnetization 

^^ = ^E(-i)'^(-<)' (5) 

i 

where (—1)* symbolizes that we sum up the contributions 
from opposite sublattices with opposite signs. This is 
the appropriate order parameter for antiferromagnetism 
when examining spin polarization occurring in bipartite 
lattices. The larger the stronger is, the antiferromag- 
netic phase. In addition to //f , the shift in the electron 
and the hole energy spectra which arises from the mag- 
netic order is quantified as AE = [E"'^^ + E^^^) /2, 
where E^'-'^ and E^^^ are the highest occupied and 
lowest unoccupied states in the ground state at half fill- 
ing, respectively. Note that in the nonmagnetic state we 
have AE = 0. We will explore how the maximum mag- 
netic moment rrimax varies with the ring width. 

The distribution of the local magnetic moments in the 
9ac '■ NzG rings for several values of N is shown in 
Fig. 1(c). The symmetry of each hexagonal ring is Cg^, 
whereas the symmetry of the magnetic moment distribu- 
tion is ICqv, i-e. the magnetic moments alter sign when 
rotated over 7r/3 rad. Therefore, it suffices to display the 
distribution of magnetic moments in sectors of 7r/3 rad, 
as done in Fig. 1(c), which combines the sectors of dif- 
ferent N. Orientation and color of triangles denotes the 
orientation of the magnetic moments, and the absolute 
value of rrii is depicted by color intensity. 
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FIG. 2; (Color online) Density of states in the 9ac '■ Nzc 
rings for A'^ taking the values 3, 4, 5, 9, 10 and 11 in the 
noninteracting system (black lines) and the interacting system 
(purple lines). 



It is evident in Fig. 1(c) that both the appearance of 
staggered magnetization and the total magnetic moment 
situated on the inner edge of the ring depend on the ring 
width. Furthermore, we observe a phase change from 
nonmagnetic order for < 3 to antiferromagnetic order 
for A'^ > 4, which is similar to previous calculations for 
zigzag hexagonal graphene dotsi^^^^ No magnetic order- 
ing for zigzag segments shorter than three unit cells is 
found because of the close proximity of the opposite sub- 
lattice imbalance on the adjacent sides of the ring inner 
edge. When this edge is short, the edge states on the 
different sides of the inner ring boundary are subject to 
strong hybridization, and therefore their energy is lifted 
from the Dirac point. Hence, spontaneous spin polar- 
ization does not occur, which is similar to the case of 
nanoribbonsj^ 

For > 4, the spatial spin symmetry is broken due 
to the e-e interaction. When the ring width decreases, 
the maximum magnetic moment, which is located near 
the middle of the zigzag edge segment increases. Fur- 
thermore, nonzero magnetization is build up on the outer 
ring edge, and it increases when the ring width decreases. 
However, as a consequence of the increasing influence of 
the outer edge with decreasing ring width, the difference 
between the distributions of the magnetic moments on 
the two edges is not large for A^ = 10 and A^ = 11. Sim- 
ilarly, the staggered magnetization increases when the 
ring width decreases. 

Figure 2 shows how the density of states (DOS) of the 




FIG. 3: (Color online) Distribution of magnetic moments in 
the 13zG : Nzg rings for A'^ ranging from 6 to 11. 



9 AC : Nzg rings (the cases depicted in Fig. 1) varies with 
A^. The density of states for the noninteracting (inter- 
acting) case is displayed by the black (purple) lines. In 
order to align the interacting and noninteracting spectra 
for easier comparison we subtracted AE for each inter- 
acting spectrum. Note that the density of states is spin 
degenerate, which is in accordance with Lieb's theorem. 
For A^ = 3, magnetic order is not present, therefore the 
energy dependence of the density of states for the inter- 
acting and noninteracting systems coincide [see Fig. 2 
(a)]. The interacting and noninteracting electron case 
exhibit a small difference in the energy dependence of 
the DOS for rings with A^ = 4 and N = 5, which is 
shown in Figs. 2(b) and 2(c). As could be inferred from 
Fig. 1(c), the magnetization along the inner ring edge 
is rather small for these values of A''. For larger A'', the 
discrepancy between the DOS's for the interacting and 
non-interacting systems becomes larger, as demonstrated 
by Figs. 2(d)-(f) for A^ = 9, 10, and 11. In all these cases, 
appreciable DOS for the noninteracting system is found 
around zero energy. Such a configuration becomes un- 
stable in the presence of e-e interactions, which results in 
the appearance of an interaction gap. 

In order to demonstrate how the shape of the outer 
boundary affects the distribution of the magnetic mo- 
ments in the ring, we show in Fig. 3 the magnetization 
in the 13^0 ■ Nzg rings. It is apparent that the shape of 
the outer edge has a large effect on the localization of the 
magnetic moments on this boundary (compare Figs. 1(c) 
and 3). It is clear that in the case of zigzag outer ring 
edge, the magnetization propagates much further into the 
ring. 

Fig. 4 displays how /i^, nimax, and AE vary with the 
length of the side of the inner ring edge expressed by the 
number A^. Along with the rings whose magnetic moment 
distributions were shown in Fig. 1(c) and Fig. 3, the case 
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FIG. 4: (Color online) (a) Staggered magnetization /i^, (b) 
maximum moment mmax and (c) energy shift AE as they 
vary with the length of the side of the inner ring edge. 



of a hexagonal graphene dot having zigzag edge, is also 
displayed in Fig. 4. Both the staggered magnetization /if 
and the energy shift AE increase with TV, i.e. with the 
size of the inner ring, except for the extremely narrow 
MzG ■ NzG rings. Interestingly, the staggered magne- 
tization in the hexagonal quantum dots does not exceed 
0.02, whereas for the IS^g ■ Nzg ring it can reach almost 
up to 0.05. The nearly twofold enlargement of the stag- 
gered magnetization could be accounted for by the double 
number of zigzag edges in the Mzg '■ Nzg ring as com- 
pared to the Nzg graphene dot. On the other hand, most 
9ac '■ Nzg rings exhibit larger staggered magnetization 
and all show larger maximum magnetic moment than the 
hexagonal graphene dot. As a matter of fact, in hexag- 
onal graphene dots the zero-energy orbitals which are 
localized along the adjacent zigzag sides of the edge are 
oriented toward each other, whereas inner zigzag edges 
in rings face away from each other. Hence, hybridization 
between the states of the two edges is larger in the for- 
mer case than in the latter case. This is why 9ac ■ Nzg 
rings turn magnetic for shorter lengths of zigzag edges 
than hexagonal dots (four versus seven, respectively). 
The decrease of rrimax with N for 13zg ■ Nzg is due to 
the more effective hybridization between the quasi-zero- 
energy states localized on the inner and outer edges of 
the ring when the ring width decreases. Hence, the elec- 
tron energy shifts from the band of zero energy states, 
and therefore magnetic ordering decays, which is mani- 
fested by a smaller rrimax in the 13zg ■ ^^zg ring than in 
the llzG dot. The shapes of the AE{N) curves shown 
in Fig. 4(c) resemble the fig{N) and mmax{N) curves in 
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FIG. 5: (Color online) (a) Contour plot of the WLDOS at 
several stages of the carving process forming the 7ac : Tza 
ring and the 1 zg hexagonal dot; the number in the upper 
left corner indicates the number of bonds cut. (b) Summed 
WLDOS in the ring, the dot, and in the whole structure as 
well as the zero energy density of states versus the number 
of bonds cut. (c) Stacked plot of the density of states; red 
depicts the densities of the stages displayed in panel (a). 



Figs. 4(a) and 4(b). 

In order to elucidate the difference between magnetic 
ordering in rings and dots, one may also analyze how 
the local density of states (LDOS) depends on the ge- 
ometry of the structure. More specifically, the spatial 
distribution of the states close to zero energy determines 
how the magnetic moments evolve when the dimensions 
of the structures varies. In order to enhance the con- 
tribution of the low-energy states, we will compute the 
weighted LDOS (WLDOS) >ii 



(6) 



Here, i indexes the lattice sites, j labels the eigenstates, 
13 is the damping coefficient chosen as l/^fP = 0.1 eV, 
whereas (pji is the value of the probability amplitude of 
the j-th state at the site i. Such defined WLDOS as- 
sumes that the contribution of the states with \Ej\ > 0.1 
eV is negligible. The plots of the WLDOS in Fig. 5(a) 
illustrate how the edge states form when the inner 7zg 
hexagonal dot is cleaved out of the outer Tag hexagonal 
dot. The inner dot is separated from the ring by severing 
the bonds one by one between the dot and the ring. The 
number of severed bonds between the dot and the ring is 
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explicitly shown in Fig. 5(a), and the dot edge is depicted 
by the blue line. The local sublattice imbalance accumu- 
lates quickly with the number of severed bonds, but no 
edge states emerge when the number of cut bonds is less 
than four. The edge states, which are depicted by red 
contours around the edge, are initially distributed evenly 
between the ring and the dot, but they extend more to 
the ring when the number of cut bonds increases. 

To explore this finding in more detail, we show in Fig. 
5(b) how the total WLDOS (full purple circles), which 
is the sum over the atomic sites in the dot (full blue cir- 
cles) and the ring (empty red circles), varies with the 
number of severed bonds. Also, the DOS at zero en- 
ergy is shown by the black triangles in Fig. 5(b). Notice 
that the variation of the WLDOS has a similar shape for 
each side of the ring's inner edge, and that the WLDOS 
displays step-like features. These steps arise because the 
imbalance between the two local sublattices, found at the 
ring and dot sides of the newly formed edge, are maxi- 
mized when the formation of each side of the rings inner 
edge is completed. The next side of the rings inner edge 
contains the opposite sublattice imbalance, and therefore 
the states on this side hybridize with the states on the 
previous side, which leads to a decrease of WLDOSJ^ 
Note that after the first edge has been cut the ring and 
the dot WLDOSs start deviating from each other more 
strongly. This is because the hybridization in the dot 
is stronger, as the edge states on adjacent segments hy- 
bridize inward and towards each other. In the ring part 
the edge states face away from each other and hybridize 
radially outward, hence the hybridization is weaker. This 
is why the WLDOS in the former case experiences a de- 
cline with the beginning of each new edge segment, while 
in the latter case the WLDOS keeps growing. The grad- 
ual increase of WLDOS for both cases near the end of 
each segment is related to the accumulation of the local 
sublattice imbalances. This pattern reappears with each 
new zigzag segment, with the exception of the last bond, 
which after being cut results in the separation of the two 
structures. By the end, the WLDOS in the ring is much 
larger than WLDOS in the dot, which accounts for the 
fact that the rings exhibit a larger maximum magnetic 
moment and staggered magnetization than the dots. Fig- 
ure 5(c) shows a stacked plot of DOSs for each resulting 
structure. Plots are stacked from the bottom up, with 
each subsequent line corresponding to a structure with 
one more bond cut. DOSs for structures depicted in Fig. 
5(a) are colored red. It shows that only features near zero 
energy evolve in a similar fashion as do the WLDOSs dur- 
ing the separation of the ring and the dot. This justifies 
the damping of states higher than 0.1 eV in calculating 
the WLDOS, as they are not artifacts of the edge forming 
between the ring and the dot. 

Finally, we examine the influence of the edge de- 
formations on somewhat larger structures; namely, the 
25zG ■ and 16 ac ■ "^Ozg rings and the 20^0 dot. 

Larger structures are considered here because they can 
be deformed in a larger number of ways than smaller 




FIG. 6: (Color online) (a) Perfect edges (black region) are 
randomly perturbed (red lines) to produce a random set of 
defects, (b) Final outlook of the deformed ring, (c) Magnetic 
moments distributions in 25zg '■ 20zg ring for several values 
of fdef- The scale is the same as in Fig. 3. (d) Staggered 
magnetization of an ensemble of randomly defected structures 
as a function of defect ratio for 25 zg '■ 20 za ring (black dots), 
16ac : 20zG ring (red dots) ring, and 20zg hexagonal dot 
(blue dots). The polynomial fitting curves are added to guide 
the eye. 



structures analyzed in the rest of the paper. Defects are 
induced by randomly deforming the polygons which out- 
line the perfect structure as shown in Fig. 6 (a). The am- 
plitude of this deformation is itself a randomly selected 
number out of a specific range and the final structure is 
made up of all atoms that are enclosed by the deformed 
outline;^ which is shown in Fig. 6(b). In order to quan- 
tify the amount of defects, the defect ratio fdef is defined 
as a fraction of the total number of the defects, which is 
a sum of the missing and the surplus sites, and the num- 
ber of the sites in the original structure. The magnetic 
moment distributions in the 25^0 : 20^0 ring for a few 
values of fdef are shown in Fig. 6 (c). Also, variation of 
the staggered magnetization with the defect fraction for 
the 16ac : 20^0 and 25^0 : 20^0 rings and the 20^0 
dot is displayed in Fig. 6 (d). For the 25^^ ■ '20 zg ring 
and the 20^0 hexagon, /if decreases with defect frac- 
tion. This is expected, having in mind that the defects 
can only impair the conditions for magnetism in zigzag 
edges. On the other hand, for the 16ac ■ 20zg ring, 
small random defects are more likely to make the larger 
outer edge magnetic than to make the smaller inner edge 
nonmagnetic. This explains the initial rise in for fdef 
up to 0.02. 

In conclusion, we predict an antiferromagnetic phase in 
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hexagonal graphene rings with zigzag inner edge within 
the mean-field Hubbard model. The distribution of mag- 
netic moments is found to strongly depend on the type of 
outer edge, and larger antiferromagnetic order is found in 
rings than in hexagonal dots. Peculiar hybridization be- 
tween the states of adjacent sides of the inner ring edge is 
found to lead to an increase of magnetization of rings with 
respect to dots. Also, the staggered magnetization and 
the maximum magnetic moment are found to be strongly 
influenced by the size and the shape of the rings. For 
wide rings, the maximum magnetic moment is largest 
when both the inner and outer edges are zigzag. But, as 
a consequence of the hybridization between the states of 



the two edges, the maximum magnetic moment in a ring 
with armchair outer edge exceeds the one for the zigzag 
outer edge when the ring width decreases. The staggered 
magnetization in both the hexagonal dots and the rings 
with zigzag outer edge is found to decrease faster than in 
the rings with armchair outer edge when the number of 
the edge defects increases. 
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